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Impact of pickup ions on the shock front nonstationarity and 
energy dissipation of the heliospheric termination shock: 
Two-dimensional full particle simulations and comparison with 

Voyager 2 observations 
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Rui Wang^ 

ABSTRACT 

The transition between the supersonic solar wind and the subsonic he- 
liosheath, the termination shock (TS), was observed by Voyager 2 (V2) on 2007 
August 31-September 1 at a distance of 84 AU from the Sun. The data reveal 
multiple crossings of a complex, quasi-perpendicular supercritical shock. These 
experimental data are a starting point for a more sophisticated analysis that in¬ 
cludes computer modeling of a shock in the presence of pickup ions (PUIs). Here, 
for the hrst time we present two-dimensional (2-D) particle-in-cell (PIC) simula¬ 
tions of the TS including PUIs self-consistently. We also report the ion velocity 
distribution across the TS using the Faraday cup data from V2. A relatively 
complete plasma and magnetic held data set from V2 gives us the opportunity 
to do a full comparison between the experimental data and PIC simulation re¬ 
sults. Our results show that: (1) The nonstationarity of the shock front is mainly 
caused by the ripples and it can exists well even at a high percentage of PUIs 
(PUI%). The width of the ripples decreases with the increasing PUI% because 
the enhanced PUI foot makes the gyro-radius of rehected solar wind ions (SWIs) 
become smaller. (2) PUIs play a key role in the energy dissipation of the TS, 
and most of the incident ion dynamic energy is transferred to the thermal energy 
of PUIs. (3) The PIC simulation indicates for the upstream parameters cho¬ 
sen for V2 conditions that the density of PUIs is about 25% and the PUIs gain 
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the larger fraction (approximately 86.6%) on average of the downstream ther¬ 
mal pressure, consistent with V2 observations near the shock. The downstream 
pressure of PUIs can very in a range from about 61.9% to 96.3% in the direction 
perpendicular to the shock normal due to the impact of shock front ripples (i.e., 
a 2-D effect). (4) The simulated composite heliosheath ion distribution function 
is a superposition of a cold core formed by transmitted SWIs, the shoulders con¬ 
tributed by the hot reflected SWIs and directly transmitted PUIs, and the wings 
of the distribution dominated by the very hot reflected PUIs. It is similar to a 
previous theoretical heliosheath ion distribution function. (5) The V2 Faraday 
cups observed the cool core of the distribution, so we see only a tip of the ice¬ 
berg. For the evolution of the cool core distribution function across the TS, the 
computed results agree reasonably well with the V2 experimental results. The 
relevance of the shock front ripples to the multiple crossings of TS observed by 
V2 is also discussed in this paper. 

Subject headings: interplanetary medium — shock waves — plasmas — Sun: 
heliosphere 


1. Introduction 

The solar wind blows outward from the Sun and forms a bubble of solar material in 
the interstellar medium. Because the interstellar plasma confines the solar wind, it has to 
become subsonic before directly interacting with the interstellar plasma, and this transition 


occurs at the heliospheric termination sh ock (TS) fjPecker et ah 120051: iRichardson et ah 


2008 


: Jokipii 

200a; 

Bnrlaga et al. 

20Q§) 


Interstellar neutrals enter the heliosphere and 


are ionized by charge exchange with the solar wind ions (SWIs ) flVa.sv1innas &: Siscoe Ill976 


Mobius et al. 

1985; 

Lee & Id 

1987; 

Galeev & Sagdeev 


19881) . They are then picked up by 


the solar wind and are named pickup ions (PUIs). After these PUIs are convected to the TS, 


influenced by the pres ence of PUIs (ILiewer fc Goldsten Il993 : 


to diffusive shock acceleration ( 

Axford et al 

19771: 

Bell 

1978 

Blandford 

19781: 

Giacalone 

2005; 

Zank et al. 

2006; 

Fisk et al. 

20061: 

Guo & Giacalone 

2010h. The TS is strongly 


Gloeckler fc Geiss 1998 


Matsukivo et al.lboo? ), whic h makes the TS very different from plaiietary and interplanetary 


shock s inside the heliosphere (ILembege et ah Il2004 iBureess et ah II2005I: iGoncharov et ah 


20141 ). Of particular interest are the microstructure and the energy dissipation of the TS. 


First, the number density of PUIs at the TS is relatively high (~ 25%) flRichardson et al 


20081: IWu et al. 1120091) and thus may greatly modify the microstructure of the shock front. 
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The TS was expected to be a boundary that is stable on a timescale of several days. The 
Voyager 2 (V2) plasma experiment observed a decrease in the solar wind speed commencing 
on about 2007 June 9, which culminated in several observations of TS crossings (named 


Decker et al. 

2008; 

Richardson et al. 

2008 


Stone et ah 


1 flBurlaga et ah 12008 


200811 ■ At least two TS crossings 


(TS-1 and TS-5) occurred when there were gaps in the tele metry. Observations of the mag¬ 
netic field structure and dynamics of the TS were reported (IBurlaga et ah 11200811 . The TS 
location depends on the solar activity, and it moves inward and outward due to changes in 
the solar wind pressure during the rising or declining phase of the solar cycle. This motion is 
on a long time scale (several years; IWhang fc Burlaga Il2000l) . The underlying mechanisms of 
how the V2 spacecraft made so many crossings in such a short time (several hours), however, 
still remain unclear. It is believed that the multiple c rossings imply motion s of the TS, which 


would be caused by the shock front nonstationarity flBurlaga et ah II2008I) . 


1987 

Hada et al. 2003 

: Ghapman et al. 

2005: Matsukivo & Scholer 

2014). 

self-exited rip- 

pies f 

Winske & Guest 

1988: Savoini & Lembege Ill99 

1 

Lem 

3ege et al. 

2004: 

Burgess & Scholer 

2007 

, pre-existing waves or turbulence 

Giacalone 

12005: 

Guo & Giacalone 

2010). which 


likely cause the unexpected motions of the TS. These nonstationarities are predicted and 
observed by numerical simulations and satellite observations, respectively. The term “self¬ 
reformation” describes a process that the particles reflected by the shock ramp accumu¬ 
late ahead of the shock and form a shock foot, which then grows and becomes the new 
ramp. The new ramp starts to reflect incident particles, and the proces s repeats. Sel f- 


refor i nation of the shock front was predicted by both hybrid simulations (IHellinger et ah 


2002: 

Lembege et al. 

2009; 

Yuan et al. 

2009; 

Tiu 

2011; 

Su et al. 

2012 

and PIG Simula- 

tions 

f Lembege & Dawson 

1987; 

Hada et al. 

2003 

Scholer et al. 

2003; 

Nishimura et al. 

2003: 

Lee et al. 

2005a 

: Yang et al. 

2OO9I, 

20R 

1) for large Mach number and low Pi shocks. 


where Pi is the ratio of the thermal pressure of ions to the magnetic pressure. Even in the 


presence of PUIs, the self-reformation can still be retrieved in some conditions fiLee et al 


2005bl: I Chapman et al. 


2005: 

Oka et al. 

2011; 

Yang et al. 

2012a: 


20141) . For the heliospheric TS, the values of Pi and the Mach number are relatively low and 

to be in the supercritical regime and proba- 
; ed rip ples are usually found 


high, respectively. The TS is ge nerally believec 


bly undergoes self-reformati on flBurlaga et al. 


i n 2-D hybrid simulations f Thomas 


1989 


2008). Self- 


exci 


Burgess fc Scholer 


f Savoini fc Lembege 1994 : Lembege et al. 12009 : Yang 


are robust in three dimensional hybrid simulations flThomas 


et al. 


2007 ) and PIC simulations 
201201). Shock front ripple s 


1989 : 


Bellinger et al. Ill996l ). 


The filament instability of the shock front fou nd in high dimensional simulations can also 


contribute to the shock front nonstationarity fjSpitkovskv l2005l: iGuo fc Giacalone 12013 
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Caprioli &: Spitkovskv II2014J) . All of the simulations above do not incl ude PUIs. By using 
the 2-D Los Alamos hybrid simulation code with PUIs, iLiu et ah (120101) studied the Alfven- 
cyclotron and mirror modes excited in the near-TS heliosheath. The impact of the PUIs on 
the shock front ripples and self-reformation is not mentioned. It is expect ed that the density 


of PUIs at the TS is of the order of 20% to 30% of th e solar wind density flRichardson et ah 


20081: IWu et ah l2009l: iMatsukivo fc Scholer 11201411 . The impact of such high percentage 
of PUIs on the shock front ripples has not been studied yet. The relevance of the shock 
front ripples to the multiple TS crossings still remains unclear. A precise description of this 
influence would require at least a 2-D full particle model. 

Second, the relative high percentage of PUIs may greatly change the dissipation mech¬ 
anism of the shock front. The V2 has a working plasma instrument and the coverage was 
sufficient to identify three crossings of th e TS. The identihed TS-3 crossing revealed an a lmost 
classical perpendicular shock structure (iBurlaga et ah II2008I: iRichardson et ah Il2008[l . It is 
clearly evidenced that the TS heats the incident SWIs very little, and the average down¬ 
stream ion temperature is m^ u ch sm aller than predicted by the MHD Rankine-Hugoniot 
conditions. IRichardson et al. (1200811 concluded that most of the solar wind energy is trans¬ 
ferred to the PUIs or othe r energetic particle s that reside in the energy range not covered 
by V2 plasma instrument. IZank et al. I ( 19961) predicted that the TS dissipation mechanism 
would favor PUIs and leave the SWIs relatively cool. They concluded that PUIs may there¬ 
fore provide the primary dissipation mechanism for a perpendicular TS wit h SWIs playing a 


secondary role. Previous 1-D hybrid simulations with multiple ion species (IWu et al. 11200911 


suggest that the density of PUIs is about 25% and the PUIs gain the lager fraction (approx¬ 
imately 90%) of the downstream thermal pressure. They dehned a downstream pressure 
ratio X = /Pt, where Pf^^ and Pt represent the thermal pressure of PUIs and the 
total ions, respectively. However, it is not clear how the upstream plasma dynamic energy 
transfers to the PUIs inside the shock front and how the downstream pressure ratio x varies 
in different locations along the shock surface (i.e., a 2-D effect). In order to understand the 
energy dissipation process and the energy partition between different species of ions, here 
we develop a 2-D full particle model of the TS. 

It is also worth studying how the velocity distribution function downstream of the TS 
looks like from PIC simulations that include PUIs. The Voyager spacecraft are making 
in situ measurements along two different trajectories in the heliosheath, but unfortunately 
they were not designed to measure the PUIs directly. Launched on 2008 October 19, IBEX 
(iMcComas et ah 1120091 1 is measuring the energetic neutral atom (ENA) flux from the bound¬ 
ary region of the heliosphere. The interpretation of ENA fluxes measured at 1 AU by IBEX 
needs knowledge of the ion velocity distribution function in the inner heliosheath. ENAs are 
created by charge exchange of interstellar neutrals with hot heliosheath protons or ions. The 
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flux of ENAs therefore d e pends sensitively on the number of hot protons downstream of the 


TS. iHeerikhuisen et ah (120081) And that the ENA flux for a k— distribution is higher than 


that for a Maxwellian proton distribution with the same temperature. This is not surprising 
because the k— distribution contains many more particles in the wings of the distribution 
function than the corresponding Maxwellian distribution. Why the heliosheath proton dis¬ 
tribution function should be like a k— distribution is, however, unclear. The answer may well 
reside in the processing of upstream PUI distribution by the TS and subsequent relaxation 
of the processed distribution in the heliosheath. 


In this paper, we use a 2-D PIC code to investigate the impact of PUIs on the shock front 
microstructure, the energy dissipation, and the downstream particle velocity distribution 
function of the TS. This paper is organized as follows. The simulation model is described in 
Section 2. We describe the simulation results in Section 3, focusing on the impact of PUIs on 
the shock front nonstationarity and the particle energy partition. In Section 4, we compare 
the simulation results with V2 observations. Finally, conclusions are drawn in Section 5. 


Simulation model 


We use a 2-D electromagnetic PIC code to simulate the evolving structures of super¬ 
critical, collisionless, perpendicular shocks with PUIs. Simulations of nonstationary shocks 


have been performed by l-D PIC codes where PUIs are included (j Chapman et ah 12005 


Matsukivo fc Scholer II 2 OIII : lYang et ah Il2012al ). In this paper, we expand our simulation 
code to 2-D in order to investigate the ripples of the shock front. Here, the control equations 
of the PIC code are Maxwell and Newton-Lorentz equations only. Due to the numerical error 
built up during PIC simulations, Guass law cannot be assumed to be satisfied all the time. 
Instead of solving Possion equation, our code solves only two curls, i.e.. Ampere and Fara- 
da y equations. A rigorous c harge conservation method for the current deposits is described 


m 


Villasenor fc Buneman fll99ll ). By using the charge conservation method, an additional 


equation dp/dt = —V ■ J should be solved at each time step to provide the current density 
to the field update. This method has been commonly used in previous PIC simulation s 


(iBuneman Il993l : iNishikawa 119971: ICai et ah l2003l: ISoitkovskv 


2008 


Rekaa et a. 


The particle data are updated by the leap-frog method (jBirdsall fc Langdon 1991 ). 


20141) . 


For the 2-D simulations, we consider a Cartesian grid {x,y). The plasma box size along 
the shock normal and shock front is = 4096A^, and Ly = 512Ay, respectively, where the 
numerical grid spacing = Ay = 0.025c/Upi = 0.25c/upe. The spatial resolution is high 
enough to resolve the microstructures of the shock front even on the electron inertial scale 


flMazelle et ah 1120101) . The physical vectors, such as velocities of particles, and electric and 
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magnetic fields E and B, have components in three directions and also spatially depend on 
X and y. A shock is p roduced by using the so-called piston method flBurgess &: Scholer 


20071: iHao et ah 112014) . in which the plasma is injected continuously from one end of the 
simulation box (a; = 0, in our case), and reflected elastically at the other end {x = L^). The 
upstream plasma has a uniform density with 16 particles per species per cell. The fractions 
of the particles of different species can be changed for different purposes (e.g., for a 25% 
PUI case, the fractions of electrons, SWIs, and PUIs are 1, 0.75 and 0.25 respectively). The 
right boundary is assumed to be a perfectly conducting barrier. The pileup of the density 
and magnetic held creates a shock propagating in the —x direction. In the 2-D simulation, 
the boundary conditions of the electromagnetic helds in the y direction are periodic, and 
the particles that move out of one end of the simulation domain in the y direction will re¬ 
enter the domain from the other end. The initial distribution functions for the SWIs and 
electrons are both Maxwellian. PUIs are injected on a thin sphere in velocity space centered 
at VinA with a radius V^hpu a s in earlier 1-D works fjChapman et ah II2005I: lYang et ah Il2012al: 


Matsukivo fc Scholer 1120141 1. The upstream Alfven speed Va is equal to 1. To reduce the 


computational time, we use an unrea listic mass ratio for ions and electrons ru i/mp, = 100, 


and light speed c = 20 Va as usual fjChapman et ah l2005l: lOka et ah l201ll) . The basic 


parameters and conhguration are as follows. The ambient magnetic held Bo is in the Y 


direction as in previous work fjWin.ske fc Quest 


1988 


(3i is 0.04 as observed by Voyager 2 fjBurlaga et ah 


Liu et ah Il2010h. The solar w ind ion 
2008; Richardson et ah 2008). The 

ni, nswi, 


injected plasma is quasi-neutral, i.e., Uj = Ue and n* = nswi + npm, where Ug 
and npui are densities of the electrons, total ions, SWIs and PUIs, respectively. 


3. Simulation results on the impact of PUIs 


In this section, the impact of PUIs on the shock front nonstationarity will be analyzed 
in detail for three cases with different percentages of PUIs (PUI%): 0% (Run A), 10% 
(Run B), and 25% (Run C). Then we concentrate on an open question of how the upstream 
ion dynamic energy transfers to the thermal energy of PUIs and SWIs and the magnetic 
energy within the shock transition layer. Furthermore, the energy partition of SWIs and 
PUIs downstream of the TS will also be computed as in previous 1-D hybrid simulations 
f Wu et ah 2009). 


First, we investigate the impact of the relative percentage of PUIs on the shock front. 
Figure [T] is an overview of Runs A (PUI%=0), B (PUI%=10), and C (PUI%=25). In each 
panel, the surface indicates the magnetic held B in the 2-D simulation domain at a time 
t = A Figure [Bi shows the shock case in the absence of PUIs. In this case, the 
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shock front is characterized by s elf-excited ripples (rnarked by a red box') as in pre vious 2-D 


PIC simulations without PUIs flSavoini fc Lembege 111994 iLembege et ah 1120091) . Typical 
structures of a supercritical perpendicular shock such as the foot (“F”), ramp (“R”) and 
overshoot (“O”) are evident in this case (see the magnetic held prohle at R = 0 marked by 
the black solid curve). The upstream Alfvenic Mach number Ma is about 4.5. A similar plot 
for Run B is shown in Figure [T]d. In this case, a weak but broad PUI foot with an average 
amplitude of ~ 1.15 Bo (marked by “PUI F”) emerges ahead of the SWI foot (marked by 
“SWI F”). The PUI foot is stationary. Figure Hh shows the corresponding plot for Run 
C. In this high PUI% (=25) case, the amplitude of the PUI foot becomes higher, but the 
amplitude of the overshoot becomes lower than that in Runs A and B. 


As can be seen from Figure [H the shock front is nonstationary even in the high PUI% 
case due to the self-excited ripples. The impact of the relative percentage of PUIs on the 
rippling shock front can be illustrated from two aspects: wave features and particle behaviors. 
Figure [2] plots the corresponding power spectrum of the fluctuating magnetic field B. The 
color shading indicates the power \B{ky)\^ which is obtained by Fourier transforming the 
values of B along the Y direction (i.e., along the shock front) at a selected position X. In 
order to show the location of the shock front, the U-averaged magnetic field R/3 has been 
overlaid on the spectrum for reference (white curve). The horizonal dashed line indicates the 
upstream value of R/3. The magnetic held huctuations are enhanced greatly form the SWI 
foot to the ramp. With the increase of PUI%, the ripple excitation region in the x direction 
(between the two vertical red lines) becomes narrower. 


To understand why the scale of ripples along Y decreases with the relative percentage 
of PUIs, we perform particle diagnosis for the cases in Figure |2l Corresponding phase space 
diagrams (A — Ua., black dots) of SWIs are plotted in the downstream rest frame (i.e., the 
simulation frame) in Figure [21 The U—averaged magnetic held B is also shown for reference 
(blue curve). A fraction of the incident SWIs coming from the left-hand side is rehected at 
the shock front, and the rehected ions become a hot SWI population when convected back 
to the downstream. The other incident SWIs are directly transmitted to the downstream 
region and form the cool core of the downstream ion velocity distribution. The shock front 
ripple is associated with the rehected SWIs. In brief, the magnetic held of the PUI foot 
ahead of the shock ramp increases with the percentage of PUIs. In a high PUI% case, the 
gyro-radius of the R SWIs becomes smaller due to the enhanced local B caused by the R 
PUIs at the foot. This leads to a narrower ripple excitation region at the shock front. 


Figure IHshows stack plots of the B prohles at a hxed Y (= 0) versus time with diherent 
PUI%. The period of the shock front nonstationarity is about 1 — 2 127^, which is close to the 


self-reformation time observed in 1-D PIC simulations flHada et ah ll2003l:IChapman et ah 














2005 


Yane et al. 

2009 

). The time scale of the first TS crossing (named TS-2) observed 

is about 29 min (jRichardson et al. 

2008 

), which is equivalent to 1.7 Our PIG 


simulation give a similar period, so TS-2 is likely caused by the self-excited rippling at the 
shock front. The time intervals of other two crossings (TS-3 and TS-4) are about 3.7 hours 
and 2.8 hours (~ 10 respectively. These crossings are likely due to th e interaction of the 


TS w ith the pre-existing waves or turbulence with a larger temporal scale (IGuo fc Giacalone 


20121 1 . 


Second, we study the energy dissipation process within the shock trans i tion l ayer and 
the resulting energy partition in the downstream region of the TS. Wu et ah fj2009 ) studied 
the energy partition of PUIs and SWIs at the TS by 1-D hybrid simulations. For compar¬ 
ison with V2 observations, they defined the energy partition of PUIs in the downstream 
as a pressure ratio y = /Pt, where Pt and Pf^^ are the downstream thermal pres¬ 
sure of the total ions and PUIs, respectively. The pressure ratio y increases as the PUI 
relative density increases. In the high percentage (25%) PUI case, the ratio y is about 
90%, which is the e nergy fraction gain for PUIs inferred from the Voyager 2 observations by 


Richardson et al. (120081) . However, 1-D hybrid simulations have the fundamental limitation 


that t he downstream h eating occurs only in the two directions perpendicular to the magnetic 
field fjWu et al. 1 120091) . 2-D hybrid simulations of quasi-perpendicular shocks without PUIs 


mirror-like fluctuations fWinske fc Quest 


sistent with observations flLiu et al. 


2006 


1988 

Lu & Wans: 

2006 

2007; 

Richardson et al. 


Hao et al. 1120141) . con 


20071) . Here, we examine 


the energy partition of PUIs in the downstream of the TS by using a 2-D PIC code which 
allows the heating to take place in all three directions. 


For reference, we present the energy dissipation at a shock without PUIs. Figure [5] shows 
the phase space plots (X — Vi^, Viy and U^) of SWIs in the shock rest frame for Run A. Note 
that the simulations are performed in the downstream rest frame so that the downstream 
flow speed |Uds| = |Us|, where Vs is the shock propagation speed. The top three panels in 
Figure |5] show that the SWIs are mainly heated along the directions perpendicular to the 
magnetic field. The dynamic pressure of the SWIs along the shock normal can be calculated 
from Pd = nmV^, where n, m, and 14 are the density, mass, and bulk velocity of ions in the 
X direction. The thermal pressure of SWIs along the shock normal can be calculated from 
Pt = nkT, where k and T are the Boltzmann constant and the temperature of the particles. 
The magnetic pressure Pb is computed from i?^/(2/io). All the pressures are normalized by 
riomQV^ of the upstream region (Figure 01), where mo and Hq are the mass and upstream 
density of the ions, respectively. A reduction of the dynamic pressure of SWIs is evident at 
the shock transition. Both of Pt and Pt increase at the the shock transition, and their sum 
is negatively correlated with the dynamic pressure Pd- Obviously, the upstream dynamic 
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energy is transferred to the SWI thermal energy and the magnetic energy of in the shock 
transition. 


The impact of PUIs on the dissipation process and the energy partition is shown in 
Figure |6] with a moderate PUI% (10). In the PUI foot region (from X = 84 c/cUpj to 
about X = 89 c/cupj), the dynamic pressures of both PUIs and SWIs start to decrease. 
Simultaneously, the thermal pressure of PUIs increases due to the reflected PUIs (Figure 
Eti-)- There are no reflected SWIs in the PUI foot (see Figure thus the thermal pressure 
of SWIs is almost unchanged in this region. The magnetic pressure Pb slightly increases at 
the PUI foot. Therefore, most of the decreased dynamic energy of the incident ions (PUIs 
plus SWIs) is transferred to the PUIs and dissipated as their thermal energy in the PUI 
foot. In the SWI foot region (from X = 89 c/upi to about X = 90.25 c/upi), the dynamic 
pressure of SWIs decreases drastically due to the reflected SWIs, while the dynamic pressure 
of PUIs is almost unchanged. The thermal pressure of PUIs and the magnetic held pressure 
increase gradually. Thus, most of the energy of SWIs is dissipated as the thermal energy 
of SWIs in the SWI foot. In the narrow ramp region (from X = 90.25 c/ujpi to the about 
X = 90.75 c/ujpi), all the pressures decrease except the magnetic pressure, which shows 
a steep increase. Thus, the dynamic energy and the thermal energy of the total ions are 
transferred to the magnetic held in this region. Figure [61:; shows the thermal pressure ratio y 
in a downstream region. Curves in diherent colors indicate the y prohles obtained in diherent 
Y locations. The ratio y can vary versus Y due to the downstream turbulence as a remnant 
of the shock front ripples (a 2-D ehect). This can not be obtained in 1-D simulations because 
of d/dY = 0. The average value of the ratio y in the current 2-D simulations is about 57.3% 
(marked by a horizonal black dashed line). 


Then we do similar calculations for the high percentage PUI (25%) case, as shown in 
Figure [71 In contrast to the low percentage PUI (10%) case, the SWIs lose almost half of 
their dynamic energy at the PUI foot (from X = 83.75 c/ujpi to about X = 88 c/ujpi). As a 
consequence, the PUIs gain more thermal energy in this region. The gains of SWI thermal 
energy and the magnetic energy are lower than those in Run B. Most of the upstream dynamic 
energy is transferred to PUIs, and the pressure of PUIs becomes much higher than that of 
SWIs. In this case, the ratio y can vary in a range from about 61.9% to 96.3% along the Y 
direction due to the remnant effect of shock front ripples. The average value of the pressure 
ratio y in the downstream regi on is about 86.6%, which is co nsistent with previous 1 -D 
hybrid models fjWu et ah [[20091) and the V2 experimental data ([Richardson et ah [[2008h . 
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4. Comparison with Voyager 2 observations 


We compare the simulation results with V2 obser vations. Figure [8^ shows the observed 
B held at the TS-3 crossing fafter iBurlaga et ah II2008I) . The identihed TS-3 crossing revealed 
an almost classical perpendicular shock structure. The plasma instrument (Faraday cups) 
on V2 worked well during the TS crossings. Figure IHb shows the observed low energy ion 
data. The average speed corresponding to the channel numbers 1, 2, 3, 4, 5, 6, 7, 8, 9, and 
10 is about 60, 90, 119, 148, 216, 256, 300, 352, and 410 km/s, respectively. The bulk speed 
of the low energy ions decreases at the foot and reaches a ver y low value after the shock 
front. The estimated upstream ion gyro-period is about 17 min (IBurlaga et ah II2008I) . Thus 
the time span of the TS-3 crossing plotted in Figure [S] is about 8.8 hl'k The V2 spacecraft 
moved outwards with a speed ~ 18 km/s, and the TS-3 front moved inwards with a speed 
~ 68 km/s. Hence, the relative speed between V2 and the shock is about 86 km/s ~ 2.28 Va, 
where Va is the Alfven speed measured in the upstream of the TS at a distance of 84 AU 
from the sun. 


In order to generate time series in PIC simulations, a virtual probe (“VP”) is used which 
records the in situ electromagnetic field and plasma information as do ne in previous work for 
Cluster crossings of the Earth’s bow shock (jScholer fc Burgess 1120061) . The VP moves from 
the upstream to the downstream of the shock in Run C (PUI%=25). The relative speed 
between the VP and the shock front is about 2.2 Va- The in situ B field seen by the VP 
is shown in Figure [HI:;. The corresponding phase space plot of the ions (SWIs plus PUIs) 
is shown in Figure [Hli. Because the low energy plasma instrument on V2 is not sensitive 
enough to see all of the ions (e.g., the wings of the ion velocity distribution), V2 can not 
directly observe the PUIs. To imitate this effect, the color bar range of Figure [Hi has been 
set from 40 to 5000. It means that the number of ions lower than 40 in the phase space 
can not be seen by the VP. Figure |8]i shows that the simulated ion distribution across the 
shock is quite similar to that observed by V2 (Figure [8 ]d). If we set the color bar range from 
0 to 5000, the VP can see all of the ions during the shock crossing (Figure |8^). It implies 
that the thermal speed of the total ions should be higher than that estimated by using V2 
plasma data. To give a more clearer view, the downstream ion velocity distribution function 
is computed from the overshoot to the far downstream of the shock. 


Figure [9^ shows the normalized total ion velocity distribution (black curve) in the 
downstream in Run C. The total heliosheath ion distribution function is a superposition of 
cold transmitted SWIs (DT SWI), hot reflected SWIs (R SWI), hot directly transmitted 
PUIs (DT PUI), and a very hot PUI population (R PUI) that is reflected by the TS and 
then convected back to the downstream. The highlighted region in red is the part expected 
to be observed by V2. Figure [9)3 shows the composite heliosheath ion distribution function 
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modeled bv IZank et al. fj2010j) . Their form of the total ion distribution function is similar 


to that obtained from our PIC simulations (Figure [9^). For reference, the blue solid and 
red dashed curves illustrates a k— distribution and a Maxwellian distribution with the same 
downstream density and temperature. Both of the modeled composite distribution and 
the K—distribution have more ions in the wings of the distribution than a corresponding 
Maxwellian one. The R SWI population is not included in Zank’s model; our simulations 
show that the contribution of R SWIs is small compared with that of R PUIs. Previous ID 
PIC simul ations have also sh own that the fraction of R SWIs decreases with the PUI% at 
the shock flYang et al. Il2012ah . Figure [Qt shows the ion velocity distribution observed by the 
Faraday cup on V2 in the heliosheath near the TS. Compared with the experimental data, 
both of velocity distribution functions obtained from the PIC simulations and Zank’s models 
imply that V2 only saw the tip of the iceberg (i.e., the cool core of the total distribution). 
This cool part of the distribution can tell the bulk velocity of the plasma, the number density 
of SWIs, and roughly the temperature of SWIs. 


5. Conclusions and Discussion 

In this paper, we use a 2-D PIC code to investigate the impact of PUIs on the nonsta- 
tionarity and energy dissipation of the TS. We summarize our main hndings below: 

1. Different from previous 1-D simulations, we have shown that the shock front rip¬ 
ples can still exist even when the relative percentage of pickup ions approaches 25%. The 
excitation of ripples is associated with the reflected solar wind ions. In a high percentage 
(25%) pickup ion case, the gyro-radius of the reflected solar wind ions becomes smaller due 
to the enhanced local B caused by the reflected pickup ions at the shock foot. This leads 
to a narrower ripple excitation region at the shock front. The period of the shock front 
nonstationarity caused by the ripples is about 1 — 2 The multiple crossings of the TS 
in a short time duration (e.g., TS-2) can be explained by the shock front nonstationarity. 

2. The energy dissipation process is examined for shocks with different percentages of 
pickup ions. At a shock in the absence of pickup ions, the transformation of wave-particle 
energy shows that the dynamic energy at the shock front is transferred to the thermal 
energy of solar wind ions and the magnetic energy. In contrast, most of the dynamic energy 
is transferred to the thermal energy of pickup ions instead of solar wind ions at a shock in 
the presence of 25% pickup ions. 

3. In order to examine the energy partition in the downstream of the shock with pickup 
ions, we compute the pressure ratio y, as in previous 1-D hybrid simulations. In the 10% 
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pickup ion case, the ratio is about 57.3%. For the 25% pickup ion case, the average value 
of the ratio x reaches about 86.6% and can vary in a range from about 61.9% to 96.3% 
along the Y direction due to the remnant effect of shock front ripples (a 2-D effect). In 1-D 
simulations, the ratio y can not vary along the Y direction because of d/dY = 0. The y- 
averaged ratio x obtained in the 2-D simulation is consistent with that obtained in previous 
1-D hybrid simulations and that estimated by the Voyager 2 experimental data. 


4. We compare our 2-D PIC simulations with the magnetic field and plasma observations 
for a typical shock crossing (TS-3). The velocity distribution of low energy solar wind ions 
resulting from PIC simulations are quit similar to the Voyager 2 observed plasma data. In 
addition, the velocity distribution of the pickup ions, which can not be directly observed by 
Voyager 2, is also predicted future observations and studies. 


5. A composite heliosheath ion distribution function is obtained in our simulation. The 
core of the distribution function is formed by the cool directly transmitted solar wind ions, 
the shoulder is contributed by the hot directly transmitted pickup ions and reflected solar 
wind ions together, and the wing of the distribution is dominated by the very hot reflected 
pickup ions. The s hape of the total distribution is similar to the theoretical model made by 


Zank et al. ([2010). Compared with the experimental data, both of the velocity distribution 


functions obtained from the PIC simulations and Zank’s model imply that Voyager 2 only 
observed the tip of the iceberg (i.e., the cool core of the total distribution). The PIC 
simulation results may help interpret the IBEX data in probing the microphysics of the 
termination shock. 


By using a fixed shock profile, iBurrows et al. (j2010l) found that multiply reflected 


pickup ions (i.e., shock surfing accelerated pickup ions) could account for the termination 
shock downstream energy gains that are generally assumed to go into the pickup ions. We 
hnd that the termination shock is nonstationary and the shock front width changes with time. 


shock prohle ( 

Zank et al. 

1996: 

Lee et al. 

1996: 

Lioatov & Zank 

1999 

2003; 

Yane et al. 

i2009h 


Shapiro &: Ucer 
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Fig. 1.— Overview of the magnetic field profile at a time t = 4 0“^ resulting from different 
Runs: (a) PUI%=0, (b) PUI%=10, and (c) PUI%=25. The color shading indicates the 
strength of the magnetic field. The shock foot, ramp and overshoot are marked by “F”, “R” 
and “O”, respectively. The foot regions due to the reflected PUIs and SWIs are labeled as 
“Pin F” and “SWI F”, respectively. The dashed black reference line shows the upstream 
value Bo- The B profiles at R = 0 are shown by black solid curves. The shock front ripples 
are marked by a red box. 
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X / {duj .) 

pi 

Fig. 2.— The power spectrum of the fluctuating magnetic held B across the shock with 0%, 
10% and 25% PUIs (from top to bottom) at f = 4 The color shading indicates the 
power {\B{Ky)\^) which is obtained by Fourier transforming the values of B along the Y 
direction at a selected position X. The vertical red lines denote the interval of the rippling 
region. The prohle 5/3 (averaged along the Y direction) is superimposed. 
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Fig. 3.— Phase space plots (Vix versus X) of SWIs at t = 4 From top to bottom, the 
relative density of PUIs is 0%, 10%, and 25%, respectively. In each panel, the Y-averaged 
B held (blue curve) is shown for reference. The ripple region is between the two vertical red 
lines. 
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Fig. 4.— Stack-plots of the B profiles (at F = 0) at different times. From left to right, the 
relative density of Fills is 0%, 10%, and 25%, respectively. 
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Fig. 5.— (a) Phase space plots {X — Vix) of the SWIs at t = 7 in Run A (PUI%=0). 
The color shading shows logarithmic distributions of the particle number, (b-c) Similar plots 
but for X — Viy and X — Viz, respectively, (d) The dynamic pressure of SWIs (black), the 
thermal pressure of SWIs (red), and the magnetic pressure (green) across the shock front. 
The sum of the SWI thermal pressure and the magnetic pressure is also shown for reference 
(red dashed). 
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Fig. 6.— (a) Phase space plots of the SWIs and PUIs at t = 7 in Run B (PUI%=10). 
The color shading shows logarithmic distributions of the particle number, (b) The dynamic 
pressure of SWIs (black) and PUIs (blue), the thermal pressure of SWIs (red) and PUIs 
(magenta), and the magnetic pressure (green) across the shock front. The black dashed 
curve indicates the total ion dynamic pressure. The red dashed curve indicates the sum of 
the total ion thermal pressure and the magnetic pressure, (c) The pressure ratio y = Pt 

computed in a downstream region (marked by “DS box” in the top panel). 
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Fig. 7.— The same format as Figure 6, but for Run C (PUI%=25). 
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Fig. 8.— (a) 48-s averages of the magnetic field strength B at the TS (TS-3, after 


Burlaga et ah 1 120081) . (b) Velocity distribution observed by Faraday cups on V2. Differ¬ 
ent channels correspond to different velocity ranges. The color shading shows logarithmic 
distributions of the current or particle count number, (c) The in situ magnetic field strength 
B seen by a virtual probe (VP) across the TS profiles obtained from Run C (PUI%=25). 
(d) Corresponding logarithmic distributions of the ions seen by the VP across the TS. Here 
the particle counts range from 40 to 5000. (e) Similar plot as in panel (d) but with a full 
count sensitivity (counts range from 0 to 5000). It shows the complete distribution of the 
total ions across the TS. 
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Fig. 9.— (a) Velocity distribution of the ions (black) selected in the “DS box” (see Figured) 
of the TS. Contributions of SWIs and PUIs to the total distribution are indicated by red and 
blue curves, respectively. The highlight region (cool core) at the peak of the distribution 
is expected to be obs erved by V2. (b) Heliosheath particle distribution (black) modeled 


by IZank et ah (120101) . The Maxwellian (red) and k, (blue) distributions with the same 
temperature and particle density are also shown for reference, (c) Cool core of the ion 
distribution observed by V2 in the downstream region of TS-3. 












